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ABSTRACT. Many different methods exist for reducing data obtained when an astronomical source is 
studied with a two-channel polarimeter, such as a WoUaston prism system. This paper presents a rigorous 
method of reducing the data from raw aperture photometry, and evaluates errors both by a statistical 
treatment, and by propagating the measured sky noise from each frame. The reduction process performs a 
hypothesis test for the presence of linear polarization. The probability of there being a non-zero polarization 
is obtained, and the best method of obtaining the normalized Stokes Parameters is discussed. Point and 
interval estimates are obtained for the degree of linear polarization, which is subject to positive bias; and 
the polarization axis is found. 
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1. Introduction 

When performing optical polarimetry of astronomical objects, we wish to answer three 
distinct, but related, physical questions. 

Firstly, is the object polarized at all? Secondly, if it is, what is the best estimate of the 
polarization? And thirdly, what confidence can we give to this measure of polarization? 

In addition to these physical questions is a presentational one: in what format should 
the results be published, so that they will be of most utility to the scientific community? 

The questions of quantifying and presenting data on linear polarization have been dis- 
cussed at length by Simmons & Stewart (1985), who note that the traditional method used 
by optical astronomers, that of Serkowski (1958), does not give the best estimate of the 
true polarization under most circumstances. Using their recommendations, I present here 
a recipe for reducing polarimetric data. 



2. Paradigm 

In this paper, I will not consider the origin of the polarization of light. It may arise from 
intrinsic polarization of the source, from interaction with the interstellar medium, or within 
Earth's atmosphere. Each of these sources represents a genuine polarization, which must 
be taken into account in explaining the measured polarization values. 

Most modern optical polarimetry systems employ a two-channel system, normally a 
Wollaston prism. Such a prism splits the incoming light into two parallel beams ('channels') 
with orthogonal polarizations - it functions as a pair of co-located linear analyzers. The 
transmission axes of the analyzers can be changed either by placing a half-wave plate before 
the prism in the optical path, and rotating this, or by rotating the actual Wollaston prism. 
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Figure 1. Reference axis, R, relative to celestial co-ordinates. 



Such a system is incapable of distinguishing circularly polarized light from unpolarized 
light, and references to 'unpolarized' light in the remainder of this paper strictly refer to 
light which is not linearly polarized; it may be totally unpolarized (i.e. randomly polarized), 
or may include a circularly polarized component. 

Where a half-wave plate is used, an anticlockwise rotation x of the waveplate results in 
an anticlockwise rotation of r/ = 2^ of the transmission axes. (For the theory of Wollaston 
prisms and wave plates, see, for instance, Chapter 8 in Hecht (1987). ) 

We will suppose that Channel 1 of the detector has a transmission axis which can be 
rotated by some angle r/ anticlockwise on the celestial sphere, relative to a reference position 
TTjQ east of north. (See Figure |l|.) The transmission axes Tl, T2, of Channels 1 and 2 are 
hence at ryo + ^ and rjo + 90° + rj respectively. 

The reference angle r/o will depend on the construction of the polarizer, and will not, 
in general, be neatly due north. For mathematical convenience in the rest of this paper, 
we will take 770 to define a reference direction, 'R', in our instrumental co-ordinate system 
and relate all other angles to it. Such instrumental angles can then be mapped on to the 
Celestial Sphere by the addition of 770- 

Since the light emerging in the two beams has traversed identical paths until reaching the 
Wollaston prism, this method of polarimetry does not suffer from the systematic errors due 
to sky fluctuation which affect single-channel polarimetry (where a single beam polarimeter 
alternately samples the two orthogonal polarizations). 

The two channels will each feed some sort of photometric array, e.g. a CCD, which will 
record a photon count. Since such images are usually built up by a process of shifting the 
image position on the array and combining the results, we will refer to a composite image 
taken in one transmission axis orientation, ?/, as a mosaic. We will denote the rate of arrival 
of photons recorded in Channel 1 and Channel 2 by ni{rj) and ^2(7?) respectively. From 
these rates, we can calculate the total intensity (I) of the source, and the difference {S) 
between the two channels: 



1(77) = ni(?7) + n2(r/). 
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= ni (7/) - 722(77). (2) 
We can also define a normalized difference: 

4V) = J^y (3) 

The purpose of tliis paper is to discuss how to interpret and present such data. 



3. Curve Fitting for p 

Suppose we have a beam of light, which has a linearly polarized component of intensity /p, 
whose electric vector points at an angle (j) anticlockwise of R. Its unpolarized component 
is of intensity When such a beam enters our detector, we can use Malus' Law (Hecht, 
1987, §8.2.1) to deduce that 

1 



and 

"-2(??) = ^4 + Ip- s\v?{4> - rf), 
from which we find 

=Iu + Ip, (4) 
and, less trivially, 

5(7?)=/p.cos[2(,/.-7?)]. (5) 
The degree of linear polarization, p, is defined by 

^p 

and so we can obtain the normalized difference by substituting Equations ^, ^ and ^ into ^: 

5(77) = p. cos[2((/) — rj)]. 

Now, if observations have been made at a number of different angles, r]j, of the trans- 
mission axis, then a series of values for rjj and Sj{r]j) will be known, and p and (j) may be 
determined by fitting a sine curve to this data, weighted by errors ag {'r]j) as necessary. This 
method has been used, for example, by di Serego Alighieri et al. (1993, §2). (Their refine- 
ment of the method allowed for the correction of the Sj {r]j ) for instrumental polarization at 
each rjj, which was necessary as they were rotating the entire camera, their system having 
no half- wave plate.) 

We note that if there is any systematic bias of Channel 1 compared to Channel 2, this 
will show up as an 7/-independent (dc) term added to the sinusoidal component when Sj{r]j) 
is fitted to the data. Such bias could arise if an object appears close to the edge of the CCD 
in one channel, for example. 
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4. The Stokes Parameters 

Polarized light is normally quantified using Stokes' parameterisation. (For basic definitions 
see, for example, Clarke, in Gehrels (ed.) (1974).) Four variables are used, but one, V, is 
only applicable to circular polarization, which a system involving only half-wave plates and 
linear analyzers cannot measure. The total intensity, /, of the light is an absolute Stokes 
Parameter. The other two parameters are defined relative to some reference axis, which in 
our case will be R, the r/o direction. Thus we define: 

Q = 5(0°) = -5(90°), 

and 

[/ = 5(45°) = -5(135°). 

Normalized Stokes' Parameters are denoted by lower case letters {q,u,v), and are found 
by dividing the raw parameters by /. We note that 5 and the normalized s can be thought 
of as a Stokes Parameter like Q oi U, generalised to an arbitrary angle - and results which 
can be derived for 5 (or s) will apply to Q and U (or q and u) as special cases. 

If the Stokes Parameters are known, then the degree and angle of polarization can be 
found: 

p = \Jq^ + n2; (7) 

(/) = ^.tan~i(n/(7), (8) 

where the signs of q and u must be inspected to determine the correct quadrant for the 
inverse tangent. Note that S (ij), Q and U must be defined as above to be consistent with 
the choice of R as Reference. 

We must now distinguish between the true values of the Stokes Parameters for a source, 
and the values which we measure in the presence of noise. We will use the subscript to 
denote the underlying values, and the subscript i for individual measured values. 

In particular, consider a source which is not polarized, so qo = uq = 0, po = 0, and 0o is 
undefined. Since the qi and Ui include noise, they will not, in general, be zero, and because 
of the form of Equation |^, pi will be a definite-positive quantity. In short, pi is a biased 
estimator for pQ. 

There is no known unbiased estimator for po, and Simmons & Stewart (1985) discuss 
at length the question of which estimator should be used. They conclude that the Stokes 
Parameters themselves are more useful than p and (p in many applications, and it is recom- 
mended, therefore, that all published polarimetric data should ideally give the normalized 
Stokes Parameters, with or without evaluation and discussion of p and (p. 

Given this preference for the Stokes Parameters it appears that one should eschew the 
curve fitting method in favour of direct evaluation of the parameters, at least when we only 
have data for the usual angles r]j = 0°, 45°, 90°, 135°. In practice, observers will take several 
observations of an object at each transmission angle. This raises the question of how best 
to combine all the measured values qi,Ui to yield a single pair of 'best estimators' for go 
and Uq - a question which is dealt with by Clarke et al. (1983) 

On the basis of this prior work and set of recommendations, it is now possible to present 
a 'recipe' for reducing polarimetric data. 
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5. Handling the Raw Data 
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5.1. SKY NOISE AND SHOT NOISE 

The raw numbers which our photometric system produces will be a set of photon count 
rates nii{r]) and n2i{r]), together with their errors, (T„^-(r/) and crn2iiv)- Errors arise from 
three sources: photon shot noise; pixel-to-pixel variations in the sky value superimposed on 
the target object; and imperfect estimation of the modal sky value to subtract from the 
image (Sterken & Manfroid, 1992; NOAO, IRAF). 

The photons arrive at the detector according to a Poisson distribution. Let the total 
integration time for a mosaic taken at a given rotation of the polarizer, be r. If the detector 
requires e photons to arrive in order for one 'count' to be registered, then the total number 
of photons incident to produce the measured signal is ner. 

Under Poisson statistics, using units of 'numbers of photons', the standard deviation of 
the number of photons arriving in this time-bin is the square root of the mean number 
arriving, viz. = ^Jner. In our detector-based count rate units, therefore, the error con- 
tributed is (Tshot = \pnfeT. Provided ner > 10, (Clarke h Cooke, 1983, §19.10) then the 
shot noise will be normally distributed, to a good approximation. 

The modal value of a sky pixel, risky can be found by considering, say, the pixel values 
in an annulus of dark sky around the object in question, an annulus which contains T) 
pixels altogether. The root-mean-square deviation of these pixels' values about the mode 
can also be found, and we will label this, cisky Hence we can estimate the error on the mode, 

o-sky/V^- 

If we perform aperture-limited photometry on our target, with an aperture of area A 
in pixels, we must subtract the modal sky level, risky, which will introduce an error 

CTskysub = A.O^ysjl^. 

Each individual pixel in the aperture will be subject to a random sky fluctuation; adding 
these in quadrature for each of the A pixels, we obtain an error ciskyfiuc = \/^-0"sky 

Ultimately, the error on the measured, normalized, intensity, is the sum in quadrature 
of the three quantities, cishot, CTskysubi and dskyfiuc- If the areas of the aperture and annulus 
are comparable, then both the second and third terms will be significant; in practice, for 
long exposure times, the first (shot) noise term will be much smaller and can be neglected. 
This is important as, unlike the sky noise, the shot noise depends on the magnitude of the 
target object itself. If its contribution to the error terms is negligible, then sky-dominated 
error terms can be compared between objects of different brightness on the same frame. 

DATA CHECK 1. For each object observed in each channel of each mosaic, the photometry 
system will have produced a count rate rii with an error, an- . For each such measurement, 
calculate ^/ni/er and verify that it is much less than (7„ . . Then one can be certain that the 
noise terms are dominated by sky noise rather than shot noise. 

5.2. TESTING FOR DC BIAS 

In practice, for each target object, we will have taken a number of mosaics at each angle 
rjj. We can immediately use each pair of intensities nii,n2i to find Ii{rjj) and Si{r]j) using 
Equations || and 
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Since the errors on the two channels are independent, we can trivially find the errors on 
both Ii{T]j) and Si{r]j); the errors turn out to be identical, and are given by: 



(TU = (^S, = xh^nj + f^na,^- (9) 



DATA CHECK 2. Take the mean value of all the Si{r]j) by summing over all the values Si 
at all angles rij; and obtain an error on this mean by combining in quadrature the error on 
each Si. If the mean value of Si{r]j), averaged over all the angles rjj, is significantly greater 
than the propagated error, then there may be some DC bias. 

Check ^ uses Si{r]j) as a measure of excess intensity in Channel 1 over Channel 2, and 
relies on the fact that there are similar numbers of observations at rjj = tj and tjj = 7/+ 90° to 
average away effects due to polarization. If, as may happen in real data gathering exercises, 
there are not identical numbers of observations at rjj = r] and rjj = r] + 90° , this could show 
up as apparent 'DC bias' in a highly polarized object. In practice, however, we are unlikely 
to encounter this combination of events; testing for bias by the above method will either 
reveal a bias much greater than the error (where the cause should be obvious when the 
original sky images are examined); or a bias consistent with the random sky noise, in which 
case we can assume that there is no significant bias. 

5.3. OBTAINING THE STOKES PARAMETERS 

Once we are satisfied that our raw data are not biased, we can proceed. At this stage in our 
data reduction, we will find it convenient to divide our set of Si{r]j) values, together with 
their associated Ii{r]j) values, into the named Stokes Parameters, 

Qi = S,{ri, = 0°) = -Siiri, = 90°) 

and 

Ui = 5,(7?, = 45°) = -5,(7?, = 135°). 

In the rest of this paper, symbols such as 5, and as^, where not followed by (r?), can be 
read as denoting 'either Qi or C/j', 'either ag^ or au/, etc.. 

DATA REDUCTION STEP 3. For each pair of data nii{rij),n2i{r]j), produce the sum, li, 
and the difference, Qi or Ui as appropriate. Using Equation |^ produce the error common 
to the sum and difference, ag. or au^. Also find the normalized difference, qi or Ui. 

In practice, for a given target object, we will have taken a small number of measurements 
of Qi and Ui - say uq and uu respectively - with individual errors obtained for each 
measurement. If the errors on the individual values are not comparable, but vary widely, 
we may need to consider taking a weighted mean. 

DATA CHECK 4. For a set of measurements of [Si,aSi), take all the measured errors, as^; 
and so find the mean error (call this if phot/' o.n'd the maximum deviation of any individual 
error from <Sphot • If the maximum deviation is large compared to the actual error, consider 
whether you need to weight the data. 
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If the deviations are large, we can weight each data point, Si, by o"5'-^^; but we will not 
pursue the subject of statistical tests on weighted means here. In practice, one normally 
finds that the noise does not vary widely between measurements. 

We have already checked (see Check |l]) that the shot noise is negligible compared with 
the sky noise terms. Therefore, the main source of variation will be the sky noise. If the 
maximum deviation of the errors from fphot is small, then we can infer that the fluctuation 
in the sky pixel values is similar in all the mosaics. 

DATA REDUCTION STEP 5. In order to carry the statistical treatment further, we must 
assume that the sky noise is normally distributed. This is standard astronomical practice. 

DATA REDUCTION STEP 6. From the sample of Stokes Parameters li, Qi and Ui, obtained 
in Step 1^ find the two means, Q and tj , with their corresponding intensities Iq and lu; 
and find the standard deviations of the two samples, -i/^g and tpjj. 

5.4. PHOTOMETRIC AND STATISTICAL ERRORS 

Since modern photometric systems can estimate the sky noise on each frame, we are faced 
throughout our data reduction sequence with a choice between two methods for handling 
errors. We can propagate the errors on individual measurements through our calculations; 
or we can use the standard deviation, tps, of the set of sample values. Si. 

In this paper, I use the symbol to denote the measured (sky-dominated) error on 
Si, and a§ for the standard error on the estimated mean, S. The standard deviation of the 
population, which is the expected error on a single measurement Si, could be denoted as, 
but above I used <?phot to make its photometric derivation obvious. 

Using statistical estimators discards the data present in the photometric noise figures 
and uses only the spread in the data points to estimate the errors. We would expect the 
statistical estimator to be of similar magnitude to the photometric error in each case; and 
a cautious approach will embrace the greater of the two errors as the better error to quote 
in each case. 

Because we may be dealing with a small sample (size I's) for some Stokes Parameter, S, 
the standard deviation of the sample, will not be the best estimator of the population 
standard deviation. The best estimator is (Clarke & Cooke, 1983, §10.5, for example): 

^^stat = J^^.i^S- (10) 

In this special case of the population standard deviation, I have used the notation £^stat 
for clarity. Conventionally, s is used for the 'best estimator' standard deviation, but this 
symbol is already in use here for a general normalized Stokes Parameter, so in this paper I 
will use the variant form of sigma, for errors derived from the sample standard deviation, 
whence ?5 = fstat, and the (statistical) standard error on the mean is 

'^^ " ~ ^s 

The mean value of our Stokes Parameter, S, is the best estimate of the true value (So) 
regardless of the size of vs- Given a choice of errors between and z.^, we will cautiously 
take the greater of the two to be the 'best' error, which we shall denote &§■ 
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DATA CHECK 7. We now have two ways of estimating the noise on a single measurement 
of a Stokes Parameter: 

• ^phot the mean sky noise level obtained from our photometry system: Check |^ 
obtains its value and verifies that the noise levels do not fluctuate greatly about this mean. 

• Statistical fluctuations in the actual values of the Stokes Parameter in question are 
quantified by i^stat; obtained by applying Equation to the data from Step 

We would expect the two noise figures to be comparable, and this can be checked in our 
data. We may also consider photometry of other objects on the same frame: Check^ shows 
us that the errors are dominated by sky noise, and Csky should be comparable between objects, 
correcting for the different apertures used: 

dsky = £x/^2Ail + A/V). 

We therefore take the best error, as, on a Stokes Parameter, S, to be the greater of EpUot 
and (Sstat- 

If our data passes the above test, then we can be reasonably confident that the statistical 
tests we will outline in the next sections will not be invalidated by noise fluctuations. 



6. Testing for Polarization 

The linear polarization of light can be thought of as a vector of length pQ and phase angle 
^0 = 200- There are two independent components to the polarization. If either Qq or Uq 
is non-zero, the light is said to be polarized. Conversely, if the light is to be described as 
unpolarized, both Qq and Uq must be shown to be zero. 

The simplest way to test whether or not our target object emits polarized light is to 
test whether the measured Stokes Parameters, Q and U, are consistent with zero. If either 
parameter is inconsistent with zero, then the source can be said to be polarized. 

To proceed, we must rely on our assumption (Step |5|) that the sky-dominated noise 
causes the raw Stokes Parameters, Qi,Ui, to be distributed normally. Then we can perform 
hypothesis testing (Clarke & Cooke, 1983, Chapters 12 and 16) for the null hypotheses that 
Qo and Uq are zero. Here, noting that the number of samples is typically small (vq ~ z^jy < 
30) we face a choice: 

• Either: assume that the sky fluctuations are normally distributed with standard devi- 
ation (Sphoti and perform hypothesis testing on the standard normal distribution with the 
statistic: 

S — 5*0 

• Or: use the variation in the Si values to estimate the population standard deviation 
fstat) and perform hypothesis testing on the Student's t distribution with us — 1 degrees of 
freedom, using the statistic: 

^_ S-Sq 

^stat/v^ 

In either case, we can perform the usual statistical test to determine whether we can 
reject the null hypothesis that '5*0 = 0', at the C5. 100% confidence level. The confidence 
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intervals for retaining the null hypothesis will be symmetrical, and will be of the forms 
— zq < z < zq and —to < t < tg. 

The values of zq and can be obtained from tables, and we define Scg to be the greater 
of ^o-^phot/v^ ^o-^stat/v^- Then the more conservative hypothesis test will reject 
that null hypothesis at the Cs.100% confidence level when \S\ > Scg- 

In such a confidence test, the probability of making a 'Type I Error', i.e. of identifying 
an unpolarized target as being polarized in one polarization sense, is simply 1 — Cs- The 
probability of correctly retaining the 'unpolarized' hypothesis is Cs- 

The probability of making a 'Type II Error' (Clarke & Cooke, 1983, §12.7) (i.e. not 
identifying a polarized target as being polarized in one polarization sense) is not trivial 
to calculate. 

Now because there are two independent senses of linear polarization, we must consider 
how to combine the results of tests on the two independent Stokes Parameters. Suppose we 
have a source which has no linear polarization. We test the two Stokes Parameters, Q and 
U, for consistency with zero at confidence levels Cq and Cu respectively. The combined 
probability of correctly retaining the null hypothesis for both channels is Cq.Cu, and that 
of making the Type I Error of rejecting the null hypothesis in either or both channels is 
1 — Cq.Cu. Hence the overall confidence of the combined test is Cq. (7(7.100%. 

Since the null hypothesis is that po = and 0o is undefined, there is no preferred direction 
in the null system, and therefore the confidence test should not prefer one channel over the 
other. Hence the test must always take place with Cq = Cjj. 

Even so, the test does not treat all angles equally; the probability of a Type II Error 
depends on the orientation of the polarization of the source. Clearly if its polarization is 
closely aligned with a transmission axis, there is a low chance of a polarization consistent 
with the null hypothesis being recorded on the aligned axis, but a much higher chance of this 
happening on the perpendicular axis. As the alignment worsens, changing 4>q while keeping 
Po constant, the probabilities for retaining the null hypothesis on the two measurement axes 
approach one another. 

Consider the case where we have taken equal numbers of measurements in the two 
channels, so vq = vu = v, and where the errors on the measurements are all of order £^phof 
Hence we can calculate zq for the null hypothesis as above. Its value will be common to 
the Q and U channels, as the noise level and the number of measurements are the same in 
both channels. 

Now suppose that the source has intensity Iq and a true non-zero polarization pQ oriented 
at position angle (pQ. Then we can write Qq = IqPq cos{2(Pq) , and Uq = lopo sm{2(pQ) . To 
generate a Type II error, a false null result must be recorded on both axes. The probability 



of a false null can be calculated for specified po and (j)o. defining zi 
probability of such a Type II error is 



IqPq 



--phot 



x=zi cos(2</>o)+2;o ry=zi sin(2<^o)+2o 



]^ rx=zicos(2(t>o)+zo f 

= ^ / / 

/TT Jx=zi cos(2(/)o)— zo Jy 



y=zi sin(20o)-2o 



exp 



1 



dx dy. 



then the 



(11) 



Clearly this probability is not independent of 4>o. 



DATA REDUCTION STEP 8. Find the 90% confidence region limits, Qqq% and Uqq%, and 
inspect whether \Q\ < ^90% ^^'^ 1^1 < ^90% ■ 

• // both Stokes Parameters fall within the limits, then the target is not shown to be 
polarized at the 81 % confidence level. In this case we can try to find polarization with some 
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lower confidence, so repeat the test for Cq = Cjj = 85%. // the null hypothesis can he 
rejected in either channel, then we have a detection at the 72.25% confidence level. There 
is probably little merit in plumbing lower confidences than this. 

• //, however, polarization is detected in one or both of the Stokes Parameters at the 
starting point of 90%, test the polarized parameters to see if the polarization remains at 
higher confidences, say 95% and 97.5%. The highest confidence with which we can reject 
the null (unpolarized) hypothesis for either Stokes Parameter should be squared to give the 
confidence with which we may claim to have detected an overall polarization. 

In our hypothesis testing, we have made the a priori assumption that all targets are 
to be assumed unpolarized until proven otherwise. This is a useful question, as we must 
ask whether our data are worth processing further - and we ask it using the raw Stokes 
Parameters, without resorting to complicated formulae. To publish useful results, howev- 
er, we must produce the normalized Stokes Parameters, together with some sort of error 
estimate, and it is this matter which we will consider next. 



7. The Normalized Stokes Parameters 

Consider a general normalized Stokes Parameter for some angle, 77: 

^ _ _ uii - n2i 
li uii + n2i 

Clarke et al. (1983) point out that the signal/noise ratio obtained by calculating 
is much better than that obtained by simply taking the mean 



~s = ^ = ^^ (12) 



* 1=1 * i=\ 



since the Equation ^ involves the taking of only one ratio, where the two terms S and / 
have better signal/noise ratios than the individual Si and li which are ratioed in Equation 



13. 

We also note that errors on S and on / are not independent of one another. We can 



write: 

ni - n2 



n\ + n2 

If we propagate through the errors on the intensities, we find: 
1 



(14) 



711+722 ^ 

In order to simplify the calculation, we recall that in Check we checked that the errors 
on all the Si (and hence li) were similar. Thus the mean error on 077e rate in one channel 
is £phot/V^- Since the number of measurements made of S* is z^5, then 
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and the error formula approximates to: 



a-s = ~s£^^o,-\J{S-^ + I~^)lys- (16) 
In practice, we will be dealing with small polarizations, so S <^ I, and knowing s from 



Equation O, then Equation 16 approximates to: 



as He^ = (17) 

As we had before with <fstat and fphot) so now we have a choice of using sky photometry 
or the statistics to estimate errors. The above method gives us the photometric error on a 
normalized Stokes' Parameter as Sphot = £p\\ot/I = ag-y/T^'-, the statistical method would 
be to take the root-mean-square deviation of the measured Sj, obtained in Step ^, about 
Clarke et a/.'s (1983) best estimator value, s: 



r- 1 



li=l 



(18) 



\/t^S - 1 ' 

DATA REDUCTION STEP 9. Following the method outlined for finding s and a^, apply 



Equations 1< and 17 to the data obtained in Step o to obtain q with Gq and u with 



DATA CHECK 10. Using q and u, compute ?g and find Sstat for both normalized Stokes 
Parameters, and compare it with ffphot eac/i case. Verify also that the errors, £x, on 
the population standard deviations for the two Stokes Parameters are similar - this should 
follow from the S -independence of Equation ^^for small q and u. 

So which error should one publish as the best estimate, on our final s — cTs or ? 
Again, a conservative approach would be to take the greater of the two in each case. 

DATA REDUCTION STEP 11. Choose the more conservative error on each normalized 
Stokes Parameter, and record the results as qicrq and uzizau. Record also the best population 
standard deviations, dq and a^- 



8. The Degree of Linear Polarization 

8.1. THE DISTRIBUTION OF THE NORMALIZED STOKES PARAMETERS 

Having obtained estimated values for q and with conservative errors, these values - 
together with the reference angle 770 - can and should be published as the most conve- 
nient form of data for colleagues to work with. It is often desired, however, to express the 
polarization not in terms of q and n, but of p and (j). 

Simmons & Stewart (1985) discuss in detail the estimation of the degree of linear polar- 
ization. Their treatment assumes that the normalized Stokes Parameters have a normal 
distribution, and that the errors on q and u are similar. This latter condition is true for 
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small polarizations (see Check 10), but before we can proceed, we must test whether the 
former condition is satisfied. 

If one assumes (Step ^) that ni and n2 are normally distributed, one can construct, 
following Clarke et al. (1983), a joint distribution for s whose parameters are the underlying 
population means (nio,n2o) and standard deviations (cri,cr2) for the photon rates nu and 
n2i- The algebra gets a little messy here, so we define three parameters, a,(3,'j: 



1 

"=2 



1 

0-1 2 



+ 



1 /I 



1 + S 



l + S 



+ 



^2o 



(19) 
(20) 

(21) 



Using these three equations, we can write the probability distribution for s as: 



P{s) 



p. exp[ 



7] 



CJl.(T2.V TT.a^.{l + s)^ 



(22) 



This can be compared to the limiting case of the normal distribution whose mean sq 
and standard error ctq are obtained by propagating the underlying means (nio,n2o) and 



standard deviations (cJi, CT2) through Equations [1^ and |l5 



Puis) 



exp[ 



2.(70^ 



ctq-v 2vr 



(23) 



We can derive an expression for the ratio R{s) = P{s)/Pn{s), which should be close to 
unity if the normalized Stokes Parameter, s, is approximately normally distributed. 



DATA CHECK 12. 

• Estimate ti-Iq and n2o using Equations [i| and [1, and the data from Step |^ Estimate 
o"! ~ (T2 — as/V^, where as is obtained from Check^. 

• Use the values of s and dg obtained in Step \l\ as the best estimates of sq and uq. 

• Hence use a computer program to calculate and plot R{s) in the domain —Sag < s < 
+3(Ts. If R(s) is close to unity throughout this domain, then we may treat the normalized 
Stokes Parameters as being normally distributed. 



8.2. POINT ESTIMATION OF p 



If the data passes Checks 10 and 12, then we can follow the method of Simmons & Stew- 
art (1985). They 'normalize' the intensity-normalized Stokes Parameters, q and u, by divid- 
ing them by their common population standard deviation, a. For clarity of notation, in a 
field where one can be discussing both probability and polarization, I will recast their for- 
mulae, such that the measured degree of polarization, normalized as required, is here given 
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in the form m = p/a; and the actual (underlying) degree of polarization, also normalized, 
is a = pq/u. It follows from the definition of p (Equation ^) that 



y + 

If (7q = au = cr, then ap = a.ypi. 

Now, Simmons & Stewart (1985) consider the case of a 'single measurement' of each of q 
and n, whereas we have found our best estimate of these parameters following the method 
of Clarke et al. (1983) However, we can consider the whole process described by Clarke et 
al. (1983) as 'a measurement', and so the treatment holds when applied to our best estimate 
of the normalized Stokes Parameters, together with the error on that estimate. 

DATA REDUCTION STEP 13. Find dp, and hence a = ap/V2, by substituting our best 
estimates of q and u and their errors (Step llip ino Equation WA Hence calculate m: 



m = \J (p + ifl /a. 

The probability distribution F{m, a) of obtaining a measured value, m, for some under- 
lying value, a, is given by the Rice distribution (Simmons & Stewart, 1985; Wardle & 
Kronberg, 1974), which is cast in the current notation using the modified Bessel function, 
lo (Boas, 1983, as defined in Ch.l2, §17): 



F{m, a) = m. exp 



.lo(ma)... (m > 0) (25) 



F{m,a) = otherwise. 



Simmons & Stewart (1985) have tested various estimators ax for bias. They find that 
when a ^ 0.7, the best estimator is the 'Maximum Likelihood Estimator', omL; which max- 
imises F{m, a) with respect to a. So oml is the solution for a of: 

a.lQ{ma) — m.Ii{ma) = 0. (26) 

If m < 1.41 then the solution of this equation is ciml = 0. 

When 0.7, the best estimator is that traditionally used by radio astronomers, e.g. 
Wardle & Kronberg (1974). In this case, the best estimator, awK) is that which maximises 
F{m, a) with respect to m, being the solution for a of: 

(1 — w?').lQ{ma) + ma.Ii{ma) = 0. (27) 

If m < 1.00 then the solution of this equation is awK = 0. 

Simmons & Stewart (1985) graph m(a) for both cases, and so show that m is a mono- 
tonically increasing function of a, and that oml < OWK < \/m. But which estimator 
should one use? Under their treatment, the selection of one of these estimators over the 
other depends on the underlying value of a; they point out that there may be good a priori 
reasons to assume greater or lesser polarizations depending upon the nature of the source. 

If we do not make any such assumptions, we can use monotonicity of m and the inequality 
oml < OWK Vm, to find two limiting cases: 
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• Let mwKmin be the solution of the Wardle & Kronberg Equation (27) for m with 
a = 0.6. Hence if m < mwKmim then cml < OWK < 0.7 and the Maximum Likehhood 
estimator is certainly the most appropriate. Calculating, we find JTiwKmin = 1-0982 <^ 1.41 
and so the Maximum Likelihood estimator will in fact be zero. 

• Let miviLmax be the solution of Maximum Likelihood Equation (^) for m with a = 0.8. 
We find mMLmax = 1.5347. Hence if m > mMLmax, then 0.7 < oml < «WKi and Wardle & 
Kronberg's estimator will clearly be the most appropriate. 

Between these two extremes, we have oml ^ 0.7 ^ awK- This presents a problem, in that 
each estimator suggests that its estimate is more appropriate than that of the other esti- 
mator. If our measured value is m-wKmin < m < mMLmax; what should we take as our best 
estimate? We could take the mean of the two estimators, but this would divide the codomain 
of d(m) into three discontinuous regions; there might be some possible polarization which 
this method could never predict! It would be better, then, to interpolate between the two 
extremes, such that in the range rrayvKmin < m < mMLmaxi 

m — mwKmin . , ""^MLmax " "T- - /r,oN 

a = .omlH -awK- (28) 

"T-MLmax " "IWKmin "J-MLmax " "T-WKmin 

If we do not know, a priori, whether a source is likely to be unpolarized, polarized to less 
than 1%, or with a greater polarization, then a would seem to be a reasonable estimator of 
the true noise-normalized polarization, and certainly better than the biased m. 

DATA REDUCTION STEP 14. Use the above criteria to find a, and hence obtain the best 
estimate, p = a. a, of the true polarization of the target. 



8.3. A CONFIDENCE INTERVAL FOR p 

As well as a point estimate for p, we would like error bars. The Rice distribution. Equation 



25, gives the probability of obtaining some m given a, and can, therefore, be used to find 
a confidence interval for the likely values of m given a. We can define two functions, C{a) 
and U{a), which give the lower and upper confidence limits for m, with some confidence 
Cp-, integrating the Rice distribution, these will satisfy: 

fm=C{a) 

F{m,a).dm = pi (29) 



m=— oo 



and 

m=+oo 

F{m,a).dm = p2 (30) 

m=U{a) 

such that 

l-Cp=pi+p2. (31) 

Such confidence intervals are non-unique, and we need to impose an additional constraint. 
We could require that the tails outside the confidence region be equal, pi = P2, but following 
Simmons &: Stewart (1985), we shall require that the confidence interval have the smallest 
possible width, in which case our additional constraint is: 

F[l({a),a]= F[C{a),a\. (32) 
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Un a-, a 



Underlying polarization, a 

Figure 2. Confidence Intervals based on the Rice Distribution. Figure adapted from Leyslion & Bales (1997). 



Prom the form of the Rice distribution, C{a) and hl{a) will be monotonically increasing 
functions of a, as shown in Figure |2|. Given a particular underlying polarization ag, the Cp 
confidence interval (mi, 7712) can be obtained by numerically solving Equations 29 thru 32 
to yield mi = vC(ao) and m2 = U{aQ). 

Now, it can be shown (Mood et al, 1974, Ch. VIII, §4.2) that the process can also be 
inverted, i.e. if we have obtained some measured value mo, then solving for mo = ^(ai) = 
C{a2) will yield a confidence interval (oi, 02), such that the confidence of a lying within this 
interval is Cp. 

Since the contours for U{a) and C{a) cut the m-axis at non-zero values of m, we must 
distinguish three cases, depending on whether or not mo lies above one or both of the 
intercepts. The values of C{0) and U{0) depend only on the confidence interval chosen; 
substituting a = into Equations 59^ thru B2 results in the pair of equations 



Cn 



exp 



£(0)^ 



and 



exp 



(33) 



C{0).exp 



:Z^(0).exp 



(34) 



A numerical solution of this pair of equations can be found for any given confidence 
interval, C^; we find that, in 67% {la) interval, £(0) = 0.4438, Z^(0) = 1.6968, while in a 
95% (2cr) interval, £(0) = 0.1094, U{Q) = 2.5048. Hence, knowing mo, and having chosen 
our desired confidence level, we can determine the interval (ai, (12) by the following criteria: 
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• mo > U{0) 

There are non-zero solutions for both U{ai) and £(02). 

• £(0) < mo < W(0) 

In this case, ai = 0, and we must solve mo = C{a2). 

• mo < £(0) 

Here, ai = 02 = 0. 

Simmons & Stewart (1985) note that the third case is formally a confidence interval of 
zero width, and suggest that this is counter-intuitive; and they go on to suggest an ad hoc 
method of obtaining a non-zero interval. However, it is perfectly reasonable to find a finite 
probability that the degree of polarization is identically zero: the source may, after all, be 
unpolarized. This can be used as the basis of estimating the probability that there is a 
non-zero underlying polarization, as will be shown in the next section. 

DATA REDUCTION STEP 15. Knowing m from Step [TJ, find the limits (01,02) appro- 
priate to confidence intervals of 67% and 95%. Hence, multiplying by a, find the confidence 
intervals on the estimated degree of polarization. The 67%o limits may he quoted as the 'error' 
on the best estimate. 



8.4. THE PROBABILITY OF THERE BEING POLARIZATION 



Consider the contour m = lA{a) on Figure ^. As defined by Equation ^ and the inversion of 
Mood et al. (1974), it divides the domain into two regions, such that there is a probability p2 
of the underlying polarization being greater than a = lA~^{mQ). There is clearly a limiting 
case where the contour cuts the m-axis at mo, hence dividing the domain into the polarized 
region a > with probability pp, and the unpolarized region with probability 1 — pp. 



Now we may substitute the Rice Distribution, Equation |2^, into Equation 3C and evaluate 
it analytically for the limiting case, a = 0: 



pp = l- exp(-mo^/2). 



(35) 



Equation ^ hence yields the probability that a measured source actually has an underlying 
polarization. 



DATA REDUCTION STEP 16. Substitute m from Step |IJ into Equation Hence quote 
the probability that the observed source is truly polarized. 



9. The Polarization Axis 

It remains to determine the axis of polarization, for which an unbiased estimate is given by 
Equation |8|. Once again, we have a choice of using the statistical or photometric errors — 
and, indeed, a choice of raw or normalized Stokes Parameters. Since 

20 = 61 = tan"^(n/g) = tan"^(r), (36) 

our first problem is to obtain the best figure for r = u/q. 

Now, as we saw in our discussion of the best normalized Stokes Parameter, it is better to 
ratio a pair of means than to take the mean of a set of ratios. We could take r = U/Q, but for 
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a very small sample, there is the danger that the mean intensity of the Q observations will 
differ from that of the U values. Therefore, we should use the normalized Stokes Parameters, 
and the least error prone estimate of the required ratio will he f = u/q, yielding (p. 
Knowing the errors on q and u, we can find the propagated error in f: 




(37) 



given the non-linear nature of the tan function, the error on (j) should be found by separately 



calculating a+ = ^ tan (f + cj' 



and o"- 



2 tan' 



Careful attention must 



be paid in the case where the error takes the phase angle across the boundary between the 
first and fourth quadrants, as the addition of ±7r to the inverse tangent may be neccessary 
to yield a sensible error in the phase angle. 



DATA REDUCTION STEP 17. Obtain 4>, the best estimate of c 



and the propagated error 



on it, = ^(IcT+l + |c7-|), using Equations 3t and 37. Add rjQ to (f> and hence quote the 
best estimate of the polarization orientation in true celestial co-ordinates. 

For the statistical error, we note that the probability distribution of observed phase 
angles, 9 = 2(j), calculated by Vinokur (1965), and quoted in Wardle & Kronberg (1974), is: 



P{e) = exp 



sin^( 



1. 

\2tt 
where 



exp 



a2 cos^{9-eo) 



+ 



a cos{9 — Oq 



27r 



+ f[acos{e-9o)] 



sign(x) 
2^ 



exp 




sign(x) .erf (x) / Vs, 



(38) 



(39) 



and erf(x) is the error function as defined in Boas (1983), Ch.ll, §9. We do not know 
a = po/a, and will have to use our best estimate, a, as obtained from Step 14. The (7^.100% 
confidence interval on the measured angle, (^1,^2)) is given by numerically solving 



in this case we choose the symmetric interval. 



(40) 



DATA REDUCTION STEP 18. Obtain the limiting values of 4) = 9/2 for confidence inter- 
vals of 67% {la) and 95% {2a). Quote the 67% limits as = ((/>2 — </>i)/2. Choose the more 
conservative error from and a^ as the best error, a^. 



10. Comparison with Other Common Techniques 



It may be instructive to note how the process of reducing polarimetric data outlined in this 
paper compares with the methods commonly used in the existing literature. The paper by 
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Simmons & Stewart (1985) gives a thorough review of five possible point estimators for the 
degree of polarisation. One of these methods is the trivial m as an estimator of a. The other 
four methods all involve the calculation of thresholds if m < 221x then ax = 0. These 
four methods are the following: 

1. Maximum Likelihood: as defined above, oml is the value of a which maximises F{m,a) 
with respect to a. Hence qml is the solution for a of Equation 26. The limit = 1-41 
is found by a numerical method. 

2. Median: Omed fixes the distribution of possible measured values such that the actual 
measured value is the median, hence J^/Z^ F{m' ,arr,ed)-dm' = 0.5. The threshold is 

aimed = 1-18, being the solution of J^,~^'"' F{m' ,0).dm' = 0.5. 

3. Serkowski's estimator: aserk fixes the distribution of possible measured values such that 
the actual measured value is the mean, hence J^,Z^ m' .F{m' ,aserk)-dm' = m. The 

threshold is mserk = 1-25 = /™, ~q°° m'.F(m', 0).dm'. 

4. Wardle & Kronberg's method: as defined above, the estimator, awKi is that which 
maximises F{m, a) with respect to m (see Equation |2^), and 221wk = 1-00. 

Simmons & Stewart (1985) note that although widely used in the optical astronomy liter- 
ature, Serkowski's estinator is not the best for either high or low polarizations; they find 
that the Wardle & Kronberg method commonly used by radio astronomers is best when 
a ^ 0.7, i.e. when the underlying polarization is high and/or the measurement noise is very 
low. The Maximum Likelihood method, superior when a ^ 0.7 (i.e. in 'difficult' conditions 
of low polarization and/or high noise), appears to be unknown in the earlier literature. 

In this paper, I have merely provided an interpolation scheme between the point esti- 
mators which they have shown to be appropriate to the 'easy' and 'difficult' measurement 
regimes. The construction of a confidence interval to estimate the error is actually indepen- 
dent of the choice of point estimator, although (as mentioned above) I believe that Simmons 
& Stewart's (1985, §3) unwillingness to 'accept sets of zero interval as confidence intervals' 
is unfounded, since physical intuition allows for the possibility of truly unpolarised sources 
(i.e. with identically zero polarizations), and their arbitrary method of avoiding zero-width 
intervals can be dispensed with. 



11. Conclusion 

The reduction of polarimetric data can seem a daunting task to the neophyte in the field. 
In this paper, I have attempted to bring together in one place the many recommendations 
made for the reduction and presentation of polarimetry, especially those of Simmons &: 
Stewart (1985), and of Clarke et al. (1983). In addition, I have suggested that it is possible 
to develop the statistical technique used by Simmons & Stewart (1985) to obtain a simple 
probability that a measured object has non-zero underlying polarization. I have also sug- 
gested that there is a form of estimator for the overall degree of linear polarization which 
is more generally applicable than either the Maximum Likelihood or the Wardle & Kron- 
berg (1974) estimators traditionally used, and which is especially relevant in cases where 
the measured data include degrees of polarization of order 0.7 times the estimated error. 

Modern computer systems can estimate the noise on each individual mosaic of a sequence 
of images; this is useful information, and is not to be discarded in favour of a crude statistical 
analysis. A recurring theme in this paper has been the comparison of the errors estimated 
from propagating the known sky noise, and from applying sampling theory to the measured 
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intensities. Bearing this in mind, I have presented here a process for data reduction in the 
form of 18 rigorous steps and checks. The recipe might be used as the basis of an automated 
data reduction process, and I hope that it will be of particular use to the researcher - 
automated or otherwise - who is attempting polarimetry for the first time. 
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Notes 

^ In an attempt to use more consistent notation, my paper uses s for the arithmetic mean of a set of 
parameters, s for a ratio of means, and s for the best (conservative) errors on certain quantities. Clarke 
et at, however, use s for the ratio of non-normalized mean Stokes Parameters, and s(l) for the arithmetic 
mean. 

^ My paper uses r) for the instrumental angle which di Serego Alighieri et al. call 0. 

^ The Wardle & Kronberg paper reproduces Vinokur's equation (my Equation BSl) but omits the factor 
'sign(a:)' from Equation ^ on the grounds (Wardle, private communication) that the probability of x falling 
in the domain a; < is negligibly small. 



